home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Libris Britannia 4
/
science library(b).zip
/
science library(b)
/
SCIENTIF
/
3929.ZIP
/
SPECTR20.EXE
/
SPECTRUM.FOR
< prev
Wrap
Text File
|
1992-04-22
|
66KB
|
2,514 lines
* SPECTRUM.FOR V 2.0
* Routine to compute Amplitude, Auto and Cross-Spectral
* Density Functions via Fast Fourier Transforms
* David E. Hess
* Fluid Flow Group - Process Measurements Division
* Chemical Science and Technology Laboratory
* National Institute of Standards and Technology
* April 15, 1992
* References: 1. Bendat, J.S. and Piersol, A.G. 1986 Random Data,
* 2nd Ed., Wiley, New York.
* 2. Harris, F.J. 1978 On the use of windows for
* harmonic analysis with the discrete fourier
* transform, Proceedings of the IEEE, 66, 51-83.
* 3. Nuttall, A.H. 1981 Some windows with very good
* sidelobe behavior, IEEE Transactions on Acoustics,
* Speech and Signal Processing, ASSP-29, 84-91.
* 4. Stearns, S.D. 1975 Digital Signal Analysis,
* Hayden Book Co., Rochelle Park, New Jersey.
* 5. Sirivat, A. 1985 Statistical Package for the
* Analysis of Data (SPAD), private communication.
* IFMAX : maximum # of files that may be processed.
* JFMAX : maximum # of files read from the command line.
* NMAX : maximum # points per record.
* NSMAX : maximum # of filter sections.
* NUMCON : # of coefficients; order of polynomial + 1
* EPS : smallest allowable number in output data.
* Phase Adjust Code parameters:
* NPTS : # pts before current point to fit
* NSD : # points in moving average of s.d. (odd # only)
* SDLIM : quit when s.d. theta .GT. SDLIM
* This code is compiled using Microsoft Fortran V 5.1. For full
* optimization, at least 603K of RAM must be available when
* compiled. For additional information refer to documentation.
* This program reads in a data file consisting of integers
* created by an A/D sampling routine and transforms the data
* to voltages or to some other variable. This other variable
* is computed using up to a fourth order polynomial; the
* coefficients must be provided in either a sequential formatted
* (Ascii) or a sequential binary file. Alternatively, the routine
* can read REAL*4 data which has already been transformed by some
* other program. The data can then be tapered in order to
* reduce side-lobe leakage. The tapering may employ a 50% overlap
* and the user has a choice of six different window functions.
* Then, the spectra are computed; the user may choose from :
* amplitude, power or cross spectra. The user may process
* up to IFMAX files; the parameters chosen for the first file
* are then applied to each of the IFMAX files. Each record may
* contain no more than NMAX points. Note that in the case of two
* channels of data per record, each channel may consist of up to
* NMAX/2 points. The speed of the routine is limited by the large
* number of I/O operations that must take place. Reduce # of points
* per record and/or the # of records per file to increase the speed
* of the routine. Alternatively, specify in the configuration file
* a RAM disk for storage of temporary files -- this also increases
* speed. The output is stored to a .PRN file with the same name
* as the input data file in order that it may be imported into
* LOTUS and plotted. In order to avoid taking the logarithm of
* zero, any number in the output data less than EPS is set equal
* to EPS. Then, the logarithm of EPS will be a finite number.
* The order of the recursive digital filters is twice the number
* of cascaded sections desired; the number of filter sections
* must be less than or equal to NSMAX. The routine will read
* up to JFMAX filenames from the command line.
INCLUDE 'SPECTRUM.FI' ! Interface statements
IMPLICIT REAL*4 (A-H,O-Z), INTEGER*2 (I-N)
INCLUDE 'SPECTRUM.FN' ! Common statements
* Allocatable arrays may not be placed within COMMON blocks
* and may not be dummy arguments within functions or subroutines.
INTEGER*2 IFIL[ALLOCATABLE](:),GAIN(0:7)
INTEGER*2 NDATA[ALLOCATABLE,HUGE](:)
REAL*4 WIN[ALLOCATABLE,HUGE](:)
REAL*4 FR[ALLOCATABLE,HUGE](:),FI[ALLOCATABLE,HUGE](:)
REAL*4 FR2[ALLOCATABLE,HUGE](:),FI2[ALLOCATABLE,HUGE](:)
REAL*4 CONST [ALLOCATABLE](:,:),B (NUMCON)
REAL*4 CONST2[ALLOCATABLE](:,:),B2(NUMCON)
DATA TNAM,WNAM,PNAM,HNAM,SNAM,DNAM
+ /'Tukey','Welch','Parzen','Hanning',
+ 'Min Side Lobe','Max Decay'/
* Initialize.
PI=2.0*ASIN(1.0)
PI2=2.0*PI
COH=.FALSE.
NOCOH=.TRUE.
ROOT1=.TRUE.
ROOT2=.TRUE.
SUFX=.FALSE.
* Code to read data from configuration file.
INQUIRE (FILE='SPECTRUM.CFG',EXIST=EXISTS)
IF (.NOT. EXISTS) GO TO 5
* File exists. Get the info.
CALL RDCFG (ROOT1,ROOT2,SUFX,PTH,TPTH,SUFFIX)
* Code to perform command line parsing.
5 NUMARGS=NARGS()
IF (NUMARGS .EQ. 1) THEN
NUMARGS=NUMARGS-1
GO TO 10
ELSE
NUMARGS=NUMARGS-1
IF (NUMARGS .GT. JFMAX) THEN
WRITE (*,'(/1X,A/1X,I1,A/1X,A)')
+ 'The routine is currently set to recognize up to ',
+ JFMAX,' filenames on the command line and will',
+ 'proceed with these filenames only.'
NUMARGS=JFMAX
ENDIF
DO I=1,NUMARGS
CALL GETARG (I,BUFFER,STATUS)
DO J=4,1,-1
FIRST=BUFFER(J:J)
IF (ICHAR(FIRST) .GE. 97 .AND. ICHAR(FIRST) .LE. 122) THEN
IHOLD=ICHAR(FIRST)-32
FIRST=CHAR(IHOLD)
BUFFER(J:J)=FIRST
ENDIF
ENDDO
IF (ICHAR(FIRST) .LT. 65 .OR. ICHAR(FIRST) .GT. 90) THEN
WRITE (*,'(/1X,A,I1,A/1X,A,A,A/)')
+ 'Filename ',I,' on the command line began with',
+ 'the nonalphabetic character ',FIRST,'.'
STOP 'Re-enter the command line correctly'
ENDIF
IF (.NOT. ROOT1) THEN
ISL=INDEX (PTH,'\',.TRUE.)
FLNM(I)=PTH(:ISL) // BUFFER // '.DAT'
ELSE
FLNM(I)=BUFFER // '.DAT'
ENDIF
INQUIRE (FILE=FLNM(I),EXIST=EXISTS)
IF (.NOT. EXISTS) THEN
IPD=INDEX (FLNM(I),'.',.TRUE.)
WRITE (*,'(/1X,A,A/)') 'Cannot find file : ',
+ FLNM(I)(:IPD+3)
STOP 'Re-enter the command line correctly'
ENDIF
ENDDO
CONSEC=.TRUE.
ARB=.FALSE.
IFSTRT=1
IFSTP=NUMARGS
IFTOT=NUMARGS
GO TO 40
ENDIF
* This controls the values that are fed to IF.
10 WRITE (*,'(/1X,A,A\)') 'Process files (C)onsecutively',
+ ' or in an (A)rbitrary order : '
READ (*,'(A)') INP
WRITE (*,'( )')
IF (INP .EQ. 'c') INP = 'C'
IF (INP .EQ. 'a') INP = 'A'
CONSEC=(INP .EQ. 'C')
ARB=(INP .EQ. 'A')
IF (.NOT. CONSEC .AND. .NOT. ARB) GO TO 10
* These statements will allow IF to be set equal to values
* stored in the IFIL array in the loop below.
IF (ARB) THEN
20 WRITE (*,'(1X,A\)') 'Enter number of files to process : '
READ (*,*) IFTOT
WRITE (*,'( )')
IF (IFTOT .GT. IFMAX) THEN
WRITE (*,'(1X,''# files must be .LE. '',I3,''.''//)') IFMAX
GO TO 20
ENDIF
ALLOCATE (IFIL(IFTOT), STAT=IERR) ! Allocate space for IFIL
IF (IERR .NE. 0)
+ STOP 'Problem allocating storage for IFIL. Aborting ...'
IFILMIN=IFMAX
IFILMAX=1
DO I=1,IFTOT
WRITE (*,'(1X,''Enter last digits of file '',I2,'' : ''\)') I
READ (*,*) IFIL(I)
IFILMIN=MIN0(IFIL(I),IFILMIN) ! Minimum of IFIL array
IFILMAX=MAX0(IFIL(I),IFILMAX) ! Maximum of IFIL array
ENDDO
WRITE (*,'( )')
IFSTRT=1
IFSTP=IFTOT
WRITE (*,'(1X,A)') 'Files are : '
WRITE (*,'(13X,10(I2,2X))') (IFIL(I),I=1,IFTOT)
WRITE (*,'( )')
ELSE
25 WRITE (*,'(1X,A\)') 'Enter starting & ending data file # : '
READ (*,*) IFSTRT,IFSTP
IF (IFSTP .LT. IFSTRT) THEN
WRITE (*,'(1X,A)')
+ 'Ending data file # >= Starting data file #.'
GO TO 25
ENDIF
WRITE (*,'( )')
IFTOT=IFSTP-IFSTRT+1
ENDIF
* Get the first letter.
30 WRITE (*,'(1X,A\)') 'Enter first letter of data file : '
READ (*,'(A)') FIRST
IF (ICHAR(FIRST) .GE. 97 .AND. ICHAR(FIRST) .LE. 122) THEN
IHOLD=ICHAR(FIRST)-32
FIRST=CHAR(IHOLD)
ENDIF
IF (ICHAR(FIRST) .LT. 65 .OR. ICHAR(FIRST) .GT. 90) THEN
WRITE (*,'(1X,A/)') 'Enter an alphabetic character (A-Z).'
GO TO 30
ENDIF
* Assign coefficient data file names.
40 IF (.NOT. ROOT1) THEN
ISL=INDEX (PTH,'\',.TRUE.)
XNAM =PTH(:ISL) // FIRST // 'CON.DAT'
YNAM =PTH(:ISL) // FIRST // 'CON.PRN'
XNAM2=PTH(:ISL) // FIRST // 'CON2.DAT'
YNAM2=PTH(:ISL) // FIRST // 'CON2.PRN'
ELSE
XNAM=FIRST//'CON.DAT'
YNAM=FIRST//'CON.PRN'
XNAM2=FIRST//'CON2.DAT'
YNAM2=FIRST//'CON2.PRN'
ENDIF
* Transform to voltages, other quantity or read Real data ?
50 WRITE (*,'(/1X,A/1X,A\)')
+ 'Transform into (V)oltages or into (O)ther quantity',
+ 'or read (R)eal data which is already transformed : '
READ (*,'(A)') INP
IF (INP .EQ. 'v') INP = 'V'
IF (INP .EQ. 'o') INP = 'O'
IF (INP .EQ. 'r') INP = 'R'
VOLT=(INP .EQ. 'V')
OTHER=(INP .EQ. 'O')
REALDAT=(INP .EQ. 'R')
IF (.NOT. VOLT .AND. .NOT. OTHER .AND.
+ .NOT. REALDAT) GO TO 50
IF (REALDAT) GO TO 70
* Get the voltage transformation factor.
60 WRITE (*,'(/1X,A\)')
+ 'Alternate voltage transformation factor (Y/N) : '
READ (*,'(A)') INP
IF (INP .EQ. 'y') INP = 'Y'
IF (INP .EQ. 'n') INP = 'N'
ALTVOLT=(INP .EQ. 'Y')
USEGAIN=(INP .EQ. 'N')
IF (.NOT. ALTVOLT .AND. .NOT. USEGAIN) GO TO 60
IF (ALTVOLT) THEN
WRITE (*,'(/1X,A\)') 'Enter voltage transformation factor : '
READ (*,*) VOFST1
ENDIF
* Remove the mean value from the file ?
70 WRITE (*,'(/1X,A,A\)') 'Remove (M)ean value,',
+ ' (D)etrend or (N)either : '
READ (*,'(A)') INP
IF (INP .EQ. 'm') INP = 'M'
IF (INP .EQ. 'd') INP = 'D'
IF (INP .EQ. 'n') INP = 'N'
RMEAN=(INP .EQ. 'M')
DTREND=(INP .EQ. 'D')
NOOP=(INP .EQ. 'N')
IF (.NOT. RMEAN .AND. .NOT. DTREND .AND. .NOT. NOOP) GO TO 70
WRITE (*,'( )')
IF (VOLT .OR. REALDAT) GO TO 90
* Determine the type of data file for coefficients.
80 WRITE (*,'(1X,A,A\)') 'Read from (A)scii or',
+ ' (B)inary data file : '
READ (*,'(A)') INP
IF (INP .EQ. 'a') INP = 'A'
IF (INP .EQ. 'b') INP = 'B'
ASCII=(INP .EQ. 'A')
BINRY=(INP .EQ. 'B')
IF (.NOT. ASCII .AND. .NOT. BINRY) GO TO 80
WRITE (*,'( )')
* Allocate space for CONST array.
IF (CONSEC) ALLOCATE (CONST(IFSTRT:IFSTP,NUMCON), STAT=IERR)
IF (ARB) ALLOCATE (CONST(IFILMIN:IFILMAX,NUMCON), STAT=IERR)
IF (IERR .NE. 0)
+ STOP 'Problem allocating storage for CONST. Aborting ...'
* Open the formatted file containing the transformation
* constants and fill the array.
IF (ASCII) THEN
INQUIRE (FILE=YNAM,EXIST=EXISTS)
IF (.NOT. EXISTS) THEN
IPD=INDEX (YNAM,'.',.TRUE.)
WRITE (*,'(/1X,A,A/)') 'Cannot find file : ',
+ YNAM(:IPD+3)
STOP ' Program terminated unsuccessfully.'
ENDIF
OPEN (NUMI,FILE=YNAM,STATUS='OLD')
READ (NUMI,*) (B(I),I=1,5)
CLOSE (NUMI,ERR=400)
DO I=IFSTRT,IFSTP
IF (CONSEC) IF=I
IF (ARB) IF=IFIL(I)
DO J=1,5
CONST(IF,J)=B(J)
ENDDO
ENDDO
ELSE IF (BINRY) THEN
* Open the binary file containing the transformation
* constants and fill the array.
INQUIRE (FILE=XNAM,EXIST=EXISTS)
IF (.NOT. EXISTS) THEN
IPD=INDEX (XNAM,'.',.TRUE.)
WRITE (*,'(/1X,A,A/)') 'Cannot find file : ',
+ XNAM(:IPD+3)
STOP ' Program terminated unsuccessfully.'
ENDIF
OPEN (NUMI,FILE=XNAM,STATUS='OLD',ACCESS='SEQUENTIAL',
+ FORM='BINARY')
READ (NUMI) NUMSETS,NSTART
IF (
+ (CONSEC .AND. (NSTART .NE. IFSTRT .OR. NUMSETS .NE. IFTOT))
+ .OR. (ARB .AND. (NSTART .NE. IFILMIN .OR.
+ NSTART+NUMSETS-1 .NE. IFILMAX)) ) THEN
WRITE (*,'(/1X,A/1X,A/1X,A/)')
+ 'Starting or ending data file #''s stored in coefficient',
+ 'data file for channel 1 do not match values entered',
+ 'into program. This must be corrected before proceeding.'
CLOSE (NUMI)
STOP ' Program terminated unsuccessfully.'
ENDIF
READ (NUMI)
+ ((CONST(NSTART+II-1,JJ),JJ=1,NUMCON),II=1,NUMSETS)
CLOSE (NUMI,ERR=400)
ENDIF
* Filter the data ?
90 WRITE (*,'(1X,A/20X,A\)')
+ 'Filter the data using (L)owpass, (H)ighpass,',
+ '(B)andpass, Band(S)top or (N)one : '
READ (*,'(A)') FTYPE
IF (FTYPE .EQ. 'l') FTYPE = 'L'
IF (FTYPE .EQ. 'h') FTYPE = 'H'
IF (FTYPE .EQ. 'b') FTYPE = 'B'
IF (FTYPE .EQ. 's') FTYPE = 'S'
IF (FTYPE .EQ. 'n') FTYPE = 'N'
LPASS=(FTYPE .EQ. 'L')
HPASS=(FTYPE .EQ. 'H')
BPASS=(FTYPE .EQ. 'B')
BSTOP=(FTYPE .EQ. 'S')
NOFIL=(FTYPE .EQ. 'N')
IF (.NOT. LPASS .AND. .NOT. HPASS .AND. .NOT. BPASS .AND.
+ .NOT. BSTOP .AND. .NOT. NOFIL) GO TO 90
IF (NOFIL) GO TO 110
* Get the cutoff frequencies.
FC=0.0
F1=0.0
F2=0.0
IF (LPASS .OR. HPASS) THEN
WRITE (*,'(/1X,A\)') 'Enter the (-3 dB) frequency cutoff : '
READ (*,*) FC
ELSE IF (BPASS .OR. BSTOP) THEN
WRITE (*,'(/1X,A\)') 'Enter the low frequency cutoff : '
READ (*,*) F1
WRITE (*,'(/1X,A\)') 'Enter the high frequency cutoff : '
READ (*,*) F2
ENDIF
100 WRITE (*,'(/1X,A/5X,A/1X,A\)')
+ 'The order of the filter will be twice',
+ 'the number of cascaded sections.',
+ 'Enter # of filter sections to cascade : '
READ (*,*) NS
IF (NS .GT. NSMAX) THEN
WRITE (*,'(1X,A,I2)') 'Reduce the # of sections to ',NSMAX
GO TO 100
ENDIF
* Two channels of data ?
110 WRITE (*,'(/1X,A\)') 'How many channels of data (1 or 2) : '
READ (*,*) ICHANS
WRITE (*,'( )')
ONECHAN=(ICHANS .EQ. 1)
TWOCHAN=(ICHANS .EQ. 2)
IF (.NOT. ONECHAN .AND. .NOT. TWOCHAN) GO TO 110
IF (ONECHAN) GO TO 170
IF (TWOCHAN .AND. REALDAT) GO TO 130
* Transform second channel to voltages or to some other quantity ?
120 WRITE (*,'(1X,A,A\)') 'Transform 2nd channel into',
+ ' (V)oltages or into (O)ther quantity : '
READ (*,'(A)') INP
IF (INP .EQ. 'v') INP = 'V'
IF (INP .EQ. 'o') INP = 'O'
VOLT2=(INP .EQ. 'V')
OTHER2=(INP .EQ. 'O')
IF (.NOT. VOLT2 .AND. .NOT. OTHER2) GO TO 120
WRITE (*,'( )')
* Get the voltage transformation factor for channel 2.
125 WRITE (*,'(1X,A\)')
+ 'Alternate voltage transformation factor (Y/N) : '
READ (*,'(A)') INP
IF (INP .EQ. 'y') INP = 'Y'
IF (INP .EQ. 'n') INP = 'N'
ALTVOLT2=(INP .EQ. 'Y')
USEGAIN2=(INP .EQ. 'N')
IF (.NOT. ALTVOLT2 .AND. .NOT. USEGAIN2) GO TO 125
WRITE (*,'( )')
IF (ALTVOLT2) THEN
WRITE (*,'(1X,A\)') 'Enter voltage transformation factor : '
READ (*,*) VOFST2
WRITE (*,'( )')
ENDIF
* Remove the mean value from the second channel of data ?
130 WRITE (*,'(1X,A,A\)') '2nd channel: remove (M)ean value,',
+ ' (D)etrend or (N)either : '
READ (*,'(A)') INP
IF (INP .EQ. 'm') INP = 'M'
IF (INP .EQ. 'd') INP = 'D'
IF (INP .EQ. 'n') INP = 'N'
RMEAN2=(INP .EQ. 'M')
DTREND2=(INP .EQ. 'D')
NOOP2=(INP .EQ. 'N')
IF (.NOT. RMEAN2 .AND. .NOT. DTREND2 .AND.
+ .NOT. NOOP2) GO TO 130
WRITE (*,'( )')
IF (VOLT2 .OR. REALDAT) GO TO 150
* Determine type of data file for coefficients for second channel.
140 WRITE (*,'(1X,A,A\)') '2nd channel: read from (A)scii or',
+ ' (B)inary data file : '
READ (*,'(A)') INP
IF (INP .EQ. 'a') INP = 'A'
IF (INP .EQ. 'b') INP = 'B'
ASCII2=(INP .EQ. 'A')
BINRY2=(INP .EQ. 'B')
IF (.NOT. ASCII2 .AND. .NOT. BINRY2) GO TO 140
WRITE (*,'( )')
* Allocate space for CONST2 array.
IF (CONSEC) ALLOCATE (CONST2(IFSTRT:IFSTP,NUMCON), STAT=IERR)
IF (ARB) ALLOCATE (CONST2(IFILMIN:IFILMAX,NUMCON), STAT=IERR)
IF (IERR .NE. 0)
+ STOP 'Problem allocating storage for CONST2. Aborting ...'
* Open the formatted file containing the transformation
* constants for the second channel and fill the array.
IF (ASCII2) THEN
INQUIRE (FILE=YNAM2,EXIST=EXISTS)
IF (.NOT. EXISTS) THEN
IPD=INDEX (YNAM2,'.',.TRUE.)
WRITE (*,'(/1X,A,A/)') 'Cannot find file : ',
+ YNAM2(:IPD+3)
STOP ' Program terminated unsuccessfully.'
ENDIF
OPEN (NUMI,FILE=YNAM2,STATUS='OLD')
READ (NUMI,*) (B2(I),I=1,5)
CLOSE (NUMI,ERR=400)
DO I=IFSTRT,IFSTP
IF (CONSEC) IF=I
IF (ARB) IF=IFIL(I)
DO J=1,5
CONST2(IF,J)=B2(J)
ENDDO
ENDDO
ELSE IF (BINRY2) THEN
* Open the binary file containing the transformation
* constants for the second channel and fill the array.
INQUIRE (FILE=XNAM2,EXIST=EXISTS)
IF (.NOT. EXISTS) THEN
IPD=INDEX (XNAM2,'.',.TRUE.)
WRITE (*,'(/1X,A,A/)') 'Cannot find file : ',
+ XNAM2(:IPD+3)
STOP ' Program terminated unsuccessfully.'
ENDIF
OPEN (NUMI,FILE=XNAM2,STATUS='OLD',ACCESS='SEQUENTIAL',
+ FORM='BINARY')
READ (NUMI) NUMSETS,NSTART
IF (
+ (CONSEC .AND. (NSTART .NE. IFSTRT .OR. NUMSETS .NE. IFTOT))
+ .OR. (ARB .AND. (NSTART .NE. IFILMIN .OR.
+ NSTART+NUMSETS-1 .NE. IFILMAX)) ) THEN
WRITE (*,'(/1X,A/1X,A/1X,A/)')
+ 'Starting or ending data file #''s stored in coefficient',
+ 'data file for channel 2 do not match values entered',
+ 'into program. This must be corrected before proceeding.'
CLOSE (NUMI)
STOP ' Program terminated unsuccessfully.'
ENDIF
READ (NUMI)
+ ((CONST2(NSTART+II-1,JJ),JJ=1,NUMCON),II=1,NUMSETS)
CLOSE (NUMI,ERR=400)
ENDIF
* Filter the data for the second channel ?
150 WRITE (*,'(1X,A/20X,A\)')
+ 'Filter the data using (L)owpass, (H)ighpass,',
+ '(B)andpass, Band(S)top or (N)one : '
READ (*,'(A)') FTYP2
IF (FTYP2 .EQ. 'l') FTYP2 = 'L'
IF (FTYP2 .EQ. 'h') FTYP2 = 'H'
IF (FTYP2 .EQ. 'b') FTYP2 = 'B'
IF (FTYP2 .EQ. 's') FTYP2 = 'S'
IF (FTYP2 .EQ. 'n') FTYP2 = 'N'
LPASS2=(FTYP2 .EQ. 'L')
HPASS2=(FTYP2 .EQ. 'H')
BPASS2=(FTYP2 .EQ. 'B')
BSTOP2=(FTYP2 .EQ. 'S')
NOFIL2=(FTYP2 .EQ. 'N')
IF (.NOT. LPASS2 .AND. .NOT. HPASS2 .AND. .NOT. BPASS2 .AND.
+ .NOT. BSTOP2 .AND. .NOT. NOFIL2) GO TO 150
WRITE (*,'( )')
IF (NOFIL2) GO TO 170
* Get the cutoff frequencies.
FC2=0.0
F12=0.0
F22=0.0
IF (LPASS2 .OR. HPASS2) THEN
WRITE (*,'(1X,A\)') 'Enter the (-3 dB) frequency cutoff : '
READ (*,*) FC2
ELSE IF (BPASS2 .OR. BSTOP2) THEN
WRITE (*,'(1X,A\)') 'Enter the low frequency cutoff : '
READ (*,*) F12
WRITE (*,'(/1X,A\)') 'Enter the high frequency cutoff : '
READ (*,*) F22
ENDIF
160 WRITE (*,'(/1X,A\)') 'Enter # of filter sections to cascade : '
READ (*,*) NS2
WRITE (*,'( )')
IF (NS2 .GT. NSMAX) THEN
WRITE (*,'(1X,A,I2)') 'Reduce the # of sections to ',NSMAX
GO TO 160
ENDIF
* Determine the kind of spectrum.
170 WRITE (*,'(1X,A\)')
+ '(A)mplitude, (P)ower or (C)ross spectrum : '
READ (*,'(A)') INP
WRITE (*,'( )')
IF (INP .EQ. 'a') INP='A'
IF (INP .EQ. 'p') INP='P'
IF (INP .EQ. 'c') INP='C'
AMP=(INP .EQ. 'A')
POWER=(INP .EQ. 'P')
CROSS=(INP .EQ. 'C')
IF (.NOT. AMP .AND. .NOT. POWER .AND. .NOT. CROSS) GO TO 170
IF (CROSS .AND. ONECHAN) THEN
WRITE (*,'(1X,A,A/)') 'Cross Spectrum requires',
+ ' two channels of data.'
GO TO 170
ENDIF
* Calculate the coherence function ?
180 IF (CROSS) THEN
WRITE (*,'(1X,A/1X,A//1X,A\)')
+ 'Coherence function requires calculation',
+ 'of both cross and autospectral density functions.',
+ 'Calculate the coherence function gamma (f) (Y/N) : '
READ (*,'(A)') INP
IF (INP .EQ. 'y') INP = 'Y'
IF (INP .EQ. 'n') INP = 'N'
COH=(INP .EQ. 'Y')
NOCOH=(INP .EQ. 'N')
IF (.NOT. COH .AND. .NOT. NOCOH) GO TO 180
WRITE (*,'( )')
ENDIF
* Taper the data or not ?
190 WRITE (*,'(1X,A\)') 'Taper the data (Y/N) : '
READ (*,'(A)') INP
WRITE (*,'( )')
IF (INP .EQ. 'y') INP='Y'
IF (INP .EQ. 'n') INP='N'
TAPER=(INP .EQ. 'Y')
NOTAPER=(INP .EQ. 'N')
IF (.NOT. TAPER .AND. .NOT. NOTAPER) GO TO 190
* Choose the window function.
IF (TAPER) THEN
200 WRITE (*,'(1X,A/1X,A\)')
+ '(H)anning, (P)arzen, (T)ukey, (W)elch,',
+ 'Min (S)ide Lobe or Max (D)ecay window : '
READ (*,'(A)') INP
WRITE (*,'( )')
IF (INP .EQ. 'h') INP='H'
IF (INP .EQ. 'p') INP='P'
IF (INP .EQ. 't') INP='T'
IF (INP .EQ. 'w') INP='W'
IF (INP .EQ. 's') INP='S'
IF (INP .EQ. 'd') INP='D'
HANN =(INP .EQ. 'H')
PARZ =(INP .EQ. 'P')
TUKEY=(INP .EQ. 'T')
WELCH=(INP .EQ. 'W')
MINSL=(INP .EQ. 'S')
MXDEC=(INP .EQ. 'D')
IF (.NOT. HANN .AND. .NOT. PARZ .AND.
+ .NOT. WELCH .AND. .NOT. TUKEY .AND.
+ .NOT. MINSL .AND. .NOT. MXDEC) GO TO 200
ENDIF
* Overlap the data or not ?
IF (TAPER) THEN
205 WRITE (*,'(1X,A\)') 'Use overlap (Y/N) : '
READ (*,'(A)') INP
WRITE (*,'( )')
IF (INP .EQ. 'y') INP='Y'
IF (INP .EQ. 'n') INP='N'
OVRLP =(INP .EQ. 'Y')
NOOVRLP=(INP .EQ. 'N')
IF (.NOT. OVRLP .AND. .NOT. NOOVRLP) GO TO 205
ENDIF
* Get the time just before starting.
CALL GETTIM (IHR,IMIN,ISEC,I100TH)
* Do the calculations for each of the desired files.
DO KOUNT=IFSTRT,IFSTP
IF (CONSEC) IF=KOUNT
IF (ARB) IF=IFIL(KOUNT)
* Assign filenames.
IF (NUMARGS .GT. 0) THEN
FIN=FLNM(IF)
ELSE
IF (.NOT. ROOT1) THEN
ISL=INDEX (PTH,'\',.TRUE.)
WRITE (FIN,'(A,A1,I3.3,A4)') PTH(:ISL),FIRST,IF,'.DAT'
ELSE
WRITE (FIN,'(A1,I3.3,A4)') FIRST,IF,'.DAT'
ENDIF
INQUIRE (FILE=FIN,EXIST=EXISTS)
IF (.NOT. EXISTS) THEN
IPD=INDEX (FIN,'.',.TRUE.)
WRITE (*,'(25X,A,A,A/)') 'Cannot find file : ',
+ FIN(:IPD+3),'.'
CYCLE
ENDIF
ENDIF
IPD=INDEX (FIN,'.',.TRUE.)
IF (ROOT2) THEN
TDOUT =FIN(IPD-4:IPD-1) // '.TRF'
TDOUT2=FIN(IPD-4:IPD-1) // '.TR2'
TOUT =FIN(IPD-4:IPD-1) // '.TAP'
TOUT2 =FIN(IPD-4:IPD-1) // '.TP2'
ELSE IF (.NOT. ROOT2) THEN
ISL=INDEX (TPTH,'\',.TRUE.)
TDOUT =TPTH(:ISL) // FIN(IPD-4:IPD-1) // '.TRF'
TDOUT2=TPTH(:ISL) // FIN(IPD-4:IPD-1) // '.TR2'
TOUT =TPTH(:ISL) // FIN(IPD-4:IPD-1) // '.TAP'
TOUT2 =TPTH(:ISL) // FIN(IPD-4:IPD-1) // '.TP2'
ENDIF
IF (.NOT. SUFX) THEN
IF (POWER) POUT=FIN(:IPD-1) // 'P.PRN'
IF (AMP) POUT=FIN(:IPD-1) // 'A.PRN'
IF (CROSS) POUT=FIN(:IPD-1) // 'C.PRN'
EOUT =FIN(:IPD-1) // 'E.PRN'
ELSE
IF (POWER) POUT=FIN(:IPD-1) // 'P' // SUFFIX
IF (AMP) POUT=FIN(:IPD-1) // 'A' // SUFFIX
IF (CROSS) POUT=FIN(:IPD-1) // 'C' // SUFFIX
EOUT =FIN(:IPD-1) // 'E' // SUFFIX
ENDIF
* Put message on screen.
IPD=INDEX (FIN,'.',.TRUE.)
IF (KOUNT .EQ. IFSTRT) THEN
WRITE (*,'(/////////////////////////13X,
+ ''S P E C T R U M C A L C U L A T I O N R O U T I N E'')')
WRITE (*,'(/25X,''Processing '',A,'' now.''/)') FIN(:IPD+3)
ELSE
WRITE (*,'(/25X,''Processing '',A,'' now.''/)') FIN(:IPD+3)
ENDIF
* Open the data file and input the data.
OPEN (NUMI,FILE=FIN,STATUS='OLD',ACCESS='SEQUENTIAL',
+ FORM='BINARY')
READ (NUMI,ERR=400) ICHANS,IRSIZE,NUMREC,IDELTMS
READ (NUMI,ERR=400) (GAIN(I),I=0,7)
* Calculate # points per record and adjust record size for
* the file containing transformed data.
IF (REALDAT) THEN
N=IRSIZE/4
IRSIZ2=IRSIZE
IRSIZE=IRSIZE/2
ELSE
N=IRSIZE/2
IRSIZ2=2*IRSIZE
ENDIF
DELT=FLOAT(IDELTMS)/1000000.0
FN=FLOAT(N)
* N less than or equal to NMAX error checking.
IF (ONECHAN) NTST=NMAX
IF (TWOCHAN) NTST=NMAX/2
IF (N .GT. NMAX) THEN
WRITE (*,'(15X,A,I5,A/15X,A,I5/)')
+ 'This data file contains ',N,' data points.',
+ '# data points per record per channel must be < ',NTST
STOP ' Program terminated unsuccessfully.'
ENDIF
* Test to make sure that N is a power of 2.
ITST=NINT(ALOG10(FN)/ALOG10(2.0))
ITST2=INT(2**ITST)-N
IF (ITST2 .NE. 0) THEN
WRITE (*,'(20X,A,I5,A/20X,A/)') 'This data file contains ',
+ N,' data points.','# data points must be a power of 2.'
STOP ' Program terminated unsuccessfully.'
ENDIF
* Allocate space for various arrays.
IF (REALDAT) THEN
ALLOCATE (WIN(N), FR (N), FI (N),
+ FR2(N), FI2(N), STAT=IERR)
IF (IERR .NE. 0)
+ STOP 'Not enough storage for data. Aborting ...'
ELSE
ALLOCATE (NDATA(N), WIN(N), FR (N), FI (N),
+ FR2(N), FI2(N), STAT=IERR)
IF (IERR .NE. 0)
+ STOP 'Not enough storage for data. Aborting ...'
ENDIF
* Test to be sure data file really has one channel.
IF (ONECHAN .AND. ICHANS .NE. 1) THEN
WRITE (*,'(13X,A,A/13X,A)') 'The header of this data',
+ ' file indicates two channels.',
+ 'You have specified one channel. Please correct problem.'
STOP ' Program terminated unsuccessfully
+.'
ENDIF
* Test to be sure data file really has two channels.
IF (TWOCHAN .AND. ICHANS .NE. 2) THEN
WRITE (*,'(13X,A,A/13X,A)') 'The header of this data',
+ ' file indicates only one channel.',
+ 'You have specified two channels. Please correct problem.'
STOP ' Program terminated unsuccessfully
+.'
ENDIF
* Estimate amount of temporary space needed.
IF (REALDAT) THEN
NBYTE=IRSIZ2*NUMREC
ELSE
NBYTE=IRSIZE*NUMREC
ENDIF
IF (.NOT. REALDAT) TRFFIL=2*NBYTE
IF (REALDAT) TRFFIL=NBYTE
IF (NOTAPER) TAPFIL=0
IF (TAPER .AND. OVRLP) TAPFIL=2*TRFFIL
IF (TAPER .AND. NOOVRLP) TAPFIL=TRFFIL
TEMP_SPACE=TRFFIL+TAPFIL
TSPCMB=FLOAT(TEMP_SPACE)/1048576.0
WRITE (*,'(13X,A,I8,A,F6.2,A/)')
+ 'Temporary space required = ',TEMP_SPACE,
+ ' bytes (',TSPCMB,' Mb).'
* Perform a write test to determine if space is available
* to hold the temporary data files.
OPEN (NUMT, FILE=TDOUT, FORM='BINARY',RECL=1,
+ ACCESS='DIRECT', STATUS='UNKNOWN')
IREC=TEMP_SPACE+(HEADER_SIZE*2)-1
WRITE (NUMT, REC=IREC, ERR=390) #FF
CLOSE (NUMT, STATUS='DELETE')
* Determine the voltage transformation factor.
IF (USEGAIN) VOFST1=(20.0/4096.0)/(2**GAIN(0))
IF (USEGAIN2) VOFST2=(20.0/4096.0)/(2**GAIN(1))
* Open the output file to hold the transformed data.
OPEN (NUMO,FILE=TDOUT,STATUS='UNKNOWN',
+ ACCESS='SEQUENTIAL',FORM='BINARY',ERR=400)
IF (TWOCHAN) OPEN (NUMO2,FILE=TDOUT2,STATUS='UNKNOWN',
+ ACCESS='SEQUENTIAL',FORM='BINARY',ERR=400)
* Write file header.
IF (ONECHAN) THEN
WRITE (NUMO,ERR=400) ICHANS,IRSIZ2,NUMREC,IDELTMS
ELSE IF (TWOCHAN) THEN
WRITE (NUMO,ERR=400) ICHANS,IRSIZE,NUMREC,IDELTMS
WRITE (NUMO2,ERR=400) ICHANS,IRSIZE,NUMREC,IDELTMS
ENDIF
* Transformation procedure for one channel of data.
IF (ONECHAN) THEN
IF (OTHER) THEN
DO J=1,NUMCON
B(J)=CONST(IF,J)
ENDDO
ELSE
B=0.0
ENDIF
IF (TRANS_ONE_CHAN (NDATA,FR,FI,B) .NE. 0) STOP
+ ' Program terminated unsuccessfully.'
* Close the input data file and proceed
* to the tapering operation.
CLOSE (NUMI,STATUS='KEEP',ERR=400)
* Transformation procedure for two channels of data.
ELSE IF (TWOCHAN) THEN
IF (OTHER) THEN
DO J=1,NUMCON
B(J)=CONST(IF,J)
ENDDO
ELSE
B=0.0
ENDIF
IF (OTHER2) THEN
DO J=1,NUMCON
B2(J)=CONST2(IF,J)
ENDDO
ELSE
B2=0.0
ENDIF
IF (TRANS_TWO_CHAN (NDATA,FR,FI,FR2,FI2,WIN,B,B2)
+ .NE. 0) STOP
+ ' Program terminated unsuccessfully.'
* Close the input data file and proceed
* to the tapering operation.
CLOSE (NUMI,STATUS='KEEP',ERR=400)
ENDIF
* Messages concerning mean removal, detrending or filtering.
MSG=.FALSE.
IF (TWOCHAN) THEN
IF (RMEAN) WRITE (*,'(25X,A,A)') 'Mean value removed',
+ ' from channel 1.'
IF (DTREND) WRITE (*,'(25X,A)') 'Channel 1 detrended.'
IF (LPASS) WRITE (*,'(25X,A,A)') 'Channel 1',
+ ' was lowpass filtered.'
IF (HPASS) WRITE (*,'(25X,A,A)') 'Channel 1',
+ ' was highpass filtered.'
IF (BPASS) WRITE (*,'(25X,A,A)') 'Channel 1',
+ ' was bandpass filtered.'
IF (BSTOP) WRITE (*,'(25X,A,A)') 'Channel 1',
+ ' was bandstop filtered.'
IF (.NOT. NOFIL) WRITE (*,'(25X,A,I2,A)')
+ 'Order of filter = ',2*NS,'.'
IF (RMEAN2) WRITE (*,'(25X,A,A)') 'Mean value removed',
+ ' from channel 2.'
IF (DTREND2) WRITE (*,'(25X,A)') 'Channel 2 detrended.'
IF (LPASS2) WRITE (*,'(25X,A,A)') 'Channel 2',
+ ' was lowpass filtered.'
IF (HPASS2) WRITE (*,'(25X,A,A)') 'Channel 2',
+ ' was highpass filtered.'
IF (BPASS2) WRITE (*,'(25X,A,A)') 'Channel 2',
+ ' was bandpass filtered.'
IF (BSTOP2) WRITE (*,'(25X,A,A)') 'Channel 2',
+ ' was bandstop filtered.'
IF (.NOT. NOFIL2) WRITE (*,'(25X,A,I2,A)')
+ 'Order of filter = ',2*NS2,'.'
IF (RMEAN .OR. DTREND .OR. .NOT. NOFIL) MSG=.TRUE.
IF (RMEAN2 .OR. DTREND2 .OR. .NOT. NOFIL) MSG=.TRUE.
ELSE
IF (RMEAN) WRITE (*,'(25X,A)') 'Mean value removed.'
IF (DTREND) WRITE (*,'(25X,A)') 'File was detrended.'
IF (LPASS) WRITE (*,'(25X,A)') 'File was lowpass filtered.'
IF (HPASS) WRITE (*,'(25X,A)') 'File was highpass filtered.'
IF (BPASS) WRITE (*,'(25X,A)') 'File was bandpass filtered.'
IF (BSTOP) WRITE (*,'(25X,A)') 'File was bandstop filtered.'
IF (.NOT. NOFIL) WRITE (*,'(25X,A,I2,A)')
+ 'Order of filter = ',2*NS,'.'
IF (RMEAN .OR. DTREND .OR. .NOT. NOFIL) MSG=.TRUE.
ENDIF
* Read the file header.
REWIND (NUMO)
READ (NUMO) ICHANS,IRSIZE,NUMREC,IDELTMS
IF (TWOCHAN) THEN
REWIND (NUMO2)
READ (NUMO2) ICHANS,IRSIZE,NUMREC,IDELTMS
ENDIF
N=IRSIZE/4
N2=N/2
NREC=NUMREC ! original numrec before overlap
NREC2=2*NUMREC
K=INT4(FLOAT(N)/10.0)
* No Tapering applied.
IF (NOTAPER) THEN
IF (MSG) THEN
WRITE (*,'(/25X,A/25X,I4.4,A)') 'No tapering applied. ',
+ NUMREC, ' records will be used.'
ELSE
WRITE (*,'(25X,A/25X,I4.4,A)') 'No tapering applied. ',
+ NUMREC, ' records will be used.'
ENDIF
GO TO 300
ENDIF
* Compute the scale factor.
IF (WIN_SCALE (WIN) .NE. 0) STOP
+ ' Program terminated unsuccessfully.'
* Open the data files to hold tapered data.
OPEN (NUMT,FILE=TOUT,STATUS='UNKNOWN',ACCESS='SEQUENTIAL',
+ FORM='BINARY',ERR=400)
IF (TWOCHAN) OPEN (NUMT2,FILE=TOUT2,STATUS='UNKNOWN',
+ ACCESS='SEQUENTIAL',FORM='BINARY',ERR=400)
IF (NOOVRLP) THEN
* No overlap - use NUMREC records.
WRITE (NUMT,ERR=400) ICHANS,IRSIZE,NUMREC,IDELTMS
IF (TWOCHAN) WRITE (NUMT2,ERR=400)
+ ICHANS,IRSIZE,NUMREC,IDELTMS
ELSE
* 50% overlap - use NREC2 records.
WRITE (NUMT,ERR=400) ICHANS,IRSIZE,NREC2,IDELTMS
IF (TWOCHAN) WRITE (NUMT2,ERR=400)
+ ICHANS,IRSIZE,NREC2,IDELTMS
ENDIF
IF (MSG) WRITE (*,'( )')
IF (HANN)
+ WRITE (*,'(25X,A)') 'Taper applied = ' // HNAM // '.'
IF (WELCH)
+ WRITE (*,'(25X,A)') 'Taper applied = ' // WNAM // '.'
IF (PARZ)
+ WRITE (*,'(25X,A)') 'Taper applied = ' // PNAM // '.'
IF (TUKEY)
+ WRITE (*,'(25X,A)') 'Taper applied = ' // TNAM // '.'
IF (MINSL)
+ WRITE (*,'(25X,A)') 'Taper applied = ' // SNAM // '.'
IF (MXDEC)
+ WRITE (*,'(25X,A)') 'Taper applied = ' // DNAM // '.'
IF (OVRLP) THEN
WRITE (*,'(25X,A/25X,I4.4,A/)') '50% overlap used.',
+ NREC2,' records will be written.'
ELSE
WRITE (*,'(25X,A/25X,I4.4,A/)') 'No overlap used.',
+ NUMREC,' records will be written.'
ENDIF
* Do the tapering now.
IF (TAPER_DATA (FR,FR2,FI,FI2,WIN) .NE. 0) STOP
+ ' Program terminated unsuccessfully.'
* Throw away transformed data file and prepare tapered data file.
CLOSE (NUMO,STATUS='DELETE',ERR=400)
IF (TWOCHAN) CLOSE (NUMO2,STATUS='DELETE',ERR=400)
REWIND (NUMT)
READ (NUMT) ICHANS,IRSIZE,NUMREC,IDELTMS
IF (TWOCHAN) THEN
REWIND (NUMT2)
READ (NUMT2) ICHANS,IRSIZE,NUMREC,IDELTMS
ENDIF
* Begin spectrum computations.
300 IF (NOTAPER) THEN
INUM=NUMO
INUM2=NUMO2
ELSE IF (TAPER) THEN
INUM=NUMT
INUM2=NUMT2
ENDIF
* Test for odd number of records. This should only occur
* at this point for the one channel, no taper or one channel,
* taper and no overlap cases.
SIZTST=FLOAT(NUMREC)
REMAIN=AMOD(SIZTST,2.0)
EVENUM=(REMAIN .EQ. 0.0)
ODDNUM=(REMAIN .NE. 0.0)
* Initialization.
N=IRSIZE/4
N2=N/2
NP1=N+1
IOD=INT4(ALOG10(FLOAT(N))/ALOG10(2.0)+.5)
NRECD2=NUMREC/2
IF (ODDNUM) NRECD2=(NUMREC+1)/2
IF (OVRLP) NUMREC=NUMREC-1
XMEAN=0.0
YMEAN=0.0
DO I=1,N
FR2(I)=0.0
FI2(I)=0.0
WIN(I)=0.0
ENDDO
IF (ONECHAN) THEN
* Spectrum computations for one channel.
IF (SPEC_ONE_CHAN (FR,FI,FR2,FI2) .NE. 0) STOP
+ ' Program terminated unsuccessfully.'
* Deallocate memory in various arrays.
IF (REALDAT) THEN
DEALLOCATE (WIN, FR, FI, FR2, FI2, STAT=IERR)
IF (IERR .NE. 0)
+ STOP 'Problem attempting to deallocate. Aborting ...'
ELSE
DEALLOCATE (NDATA, WIN, FR, FI, FR2, FI2, STAT=IERR)
IF (IERR .NE. 0)
+ STOP 'Problem attempting to deallocate. Aborting ...'
ENDIF
CYCLE
ELSE IF (TWOCHAN) THEN
* Spectrum computations for two channels.
IF (SPEC_TWO_CHAN (FR,FI,FR2,FI2,WIN) .NE. 0) STOP
+ ' Program terminated unsuccessfully.'
* Deallocate memory in various arrays.
IF (REALDAT) THEN
DEALLOCATE (WIN, FR, FI, FR2, FI2, STAT=IERR)
IF (IERR .NE. 0)
+ STOP 'Problem attempting to deallocate. Aborting ...'
ELSE
DEALLOCATE (NDATA, WIN, FR, FI, FR2, FI2, STAT=IERR)
IF (IERR .NE. 0)
+ STOP 'Problem attempting to deallocate. Aborting ...'
ENDIF
ENDIF
ENDDO
* Get the time just before ending.
CALL GETTIM (IHR2,IMIN2,ISEC2,I100TH2)
WRITE (*,'(/25X,A,I2.2,A,I2.2,A,I2.2,A,I2.2)')
+ 'Starting time : ',IHR,':',IMIN,':',ISEC,'.',I100TH
WRITE (*,'(25X,A,I2.2,A,I2.2,A,I2.2,A,I2.2)')
+ 'Ending time : ',IHR2,':',IMIN2,':',ISEC2,'.',I100TH2
* Calculate elapsed time - hundredths of a second.
IF (I100TH2-I100TH .LT. 0) THEN
I100TH2=I100TH2+100
I100H=I100TH2-I100TH
ISEC2=ISEC2-1
ELSE
I100H=I100TH2-I100TH
ENDIF
* Seconds.
IF (ISEC2-ISEC .LT. 0) THEN
ISEC2=ISEC2+60
ISECH=ISEC2-ISEC
IMIN2=IMIN2-1
ELSE
ISECH=ISEC2-ISEC
ENDIF
* Minutes.
IF (IMIN2-IMIN .LT. 0) THEN
IMIN2=IMIN2+60
IMINH=IMIN2-IMIN
IHR2=IHR2-1
ELSE
IMINH=IMIN2-IMIN
ENDIF
* Hours.
IF (IHR2-IHR .LT. 0) THEN
IHR2=IHR2+24
IHRH=IHR2-IHR
ELSE
IHRH=IHR2-IHR
ENDIF
* Display elapsed time.
WRITE (*,'(25X,A,I2.2,A,I2.2,A,I2.2,A,I2.2)')
+ 'Elapsed time : ',IHRH,':',IMINH,':',ISECH,'.',I100H
WRITE (*,'( )')
STOP ' Program terminated successfully.'
* I/O Error checking.
390 WRITE (*,'(/17X,A/)')
+ 'Not enough space available for temporary files.'
STOP ' Program terminated unsuccessfully.'
400 WRITE (*,'(/13X,A,A/)') 'Run-time I/O error : ',
+ 'No such device or out of space !'
STOP ' Program terminated unsuccessfully.'
END
FUNCTION TRANS_ONE_CHAN (NDATA,FR,FI,B)
IMPLICIT REAL*4 (A-H,O-Z), INTEGER*2 (I-N)
INCLUDE 'SPECTRUM.FN' ! Common statements
INTEGER*2 NDATA[HUGE](*)
REAL*4 FR[HUGE](*)
REAL*4 FI[HUGE](*)
REAL*4 B(*)
* Initialization.
TRANS_ONE_CHAN=0
SUM1UT=0.0
SUM2UT=0.0
* Perform mean and RMS value calculations for each record
* in this data file.
DO IB=1,NUMREC
* Fill the array with the current record.
IF (REALDAT) THEN
READ (NUMI) (FR(I),I=1,N)
DO I=1,N
FI(I)=FR(I)*FR(I)
ENDDO
ELSE
READ (NUMI) (NDATA(I),I=1,N)
ENDIF
* Transform the data into voltages or velocities.
* Calculate squared values for RMS calculations.
IF (REALDAT) GO TO 10
DO I=1,N
V=FLOAT(NDATA(I))*VOFST1
IF (VOLT) THEN
FR(I)=V
FI(I)=FR(I)*FR(I)
CYCLE
ENDIF
FR(I)=B(5)
DO J=4,1,-1
FR(I)=FR(I)*V+B(J)
ENDDO
FI(I)=FR(I)*FR(I)
ENDDO
* Determine mean and RMS values for each record.
10 SUM1U=ADDSUM(FR,N)
SUM1UT=SUM1UT+SUM1U
SUM2U=ADDSUM(FI,N)
SUM2UT=SUM2UT+SUM2U
UBAR=SUM1U/FN
URMS=SQRT(SUM2U/FN-UBAR**2)
* Display output for each record.
IF (IB .EQ. 1) THEN
WRITE (*,20) IB,UBAR,URMS
20 FORMAT (16X,'Record ',I4.4,4X,'Ave = ',G11.5,
+ 2X,'Rms = ',G11.5)
ELSE
WRITE (*,30) IB,UBAR,URMS
30 FORMAT ('+',15X,'Record ',I4.4,4X,'Ave = ',G11.5,
+ 2X,'Rms = ',G11.5)
ENDIF
* Save transformed values for this record.
IF (.NOT. NOFIL)
+ CALL FILTER (FTYPE,FC,F1,F2,NS,DELT,N,FR)
IF (RMEAN) THEN
DO I=1,N
FR(I)=FR(I)-UBAR
ENDDO
ENDIF
IF (DTREND) CALL DTRD (FR,N,DELT)
WRITE (NUMO,ERR=40) (FR(I), I=1,N)
ENDDO
* All records of this data file have been processed
* and stored. Compute mean and rms values for this
* data file and continue.
TOTAL=FLOAT(NUMREC*N)
AVEU=SUM1UT/TOTAL
RMSU=SQRT(SUM2UT/TOTAL-AVEU**2)
* Display the results for this file.
WRITE (*,'(/16X,A,1X,A,G11.5,2X,A,G11.5/)') 'File Totals : ',
+ 'Ave = ',AVEU,'Rms = ',RMSU
RETURN
40 WRITE (*,'(/13X,A,A/)') 'Run-time I/O error : ',
+ 'No such device or out of space !'
TRANS_ONE_CHAN=1
RETURN
END
FUNCTION TRANS_TWO_CHAN (NDATA,FR,FI,FR2,FI2,WIN,B,B2)
IMPLICIT REAL*4 (A-H,O-Z), INTEGER*2 (I-N)
INCLUDE 'SPECTRUM.FN' ! Common statements
INTEGER*2 NDATA[HUGE](*)
REAL*4 FR[HUGE](*)
REAL*4 FR2[HUGE](*)
REAL*4 FI[HUGE](*)
REAL*4 FI2[HUGE](*)
REAL*4 WIN[HUGE](*)
REAL*4 B(*),B2(*)
* Initialization.
TRANS_TWO_CHAN=0
SUM1UT=0.0
SUM2UT=0.0
SUM1VT=0.0
SUM2VT=0.0
N=N/2
FN=FLOAT(N)
* Perform mean and RMS value calculations for each record
* in this data file.
DO IB=1,NUMREC
* Fill the array with the current record.
IF (REALDAT) THEN
READ (NUMI) (WIN(I),I=1,2*N)
DO II=1,2*N,2
I=(II+1)/2
FR(I)=WIN(II)
FR2(I)=WIN(II+1)
ENDDO
DO I=1,N
FI(I)=FR(I)*FR(I)
FI2(I)=FR2(I)*FR2(I)
ENDDO
ELSE
READ (NUMI) (NDATA(I),I=1,2*N)
ENDIF
* Transform the data into voltages or velocities.
* Calculate squared values for RMS calculations.
IF (REALDAT) GO TO 20
DO II=1,2*N,2
I=(II+1)/2
V=FLOAT(NDATA(II))*VOFST1
V2=FLOAT(NDATA(II+1))*VOFST2
IF (VOLT) THEN
FR(I)=V
FI(I)=FR(I)*FR(I)
GO TO 10
ENDIF
FR(I)=B(5)
DO J=4,1,-1
FR(I)=FR(I)*V+B(J)
ENDDO
FI(I)=FR(I)*FR(I)
10 IF (VOLT2) THEN
FR2(I)=V2
FI2(I)=FR2(I)*FR2(I)
CYCLE
ENDIF
FR2(I)=B2(5)
DO J=4,1,-1
FR2(I)=FR2(I)*V2+B2(J)
ENDDO
FI2(I)=FR2(I)*FR2(I)
ENDDO
* Determine mean and RMS values for each record.
20 SUM1U=ADDSUM(FR,N)
SUM1UT=SUM1UT+SUM1U
SUM2U=ADDSUM(FI,N)
SUM2UT=SUM2UT+SUM2U
UBAR=SUM1U/FN
URMS=SQRT(SUM2U/FN-UBAR**2)
SUM1V=ADDSUM(FR2,N)
SUM1VT=SUM1VT+SUM1V
SUM2V=ADDSUM(FI2,N)
SUM2VT=SUM2VT+SUM2V
VBAR=SUM1V/FN
VRMS=SQRT(SUM2V/FN-VBAR**2)
* Display output for each record.
IF (IB .EQ. 1) THEN
WRITE (*,30) IB,UBAR,URMS,VBAR,VRMS
30 FORMAT (1X,'Rec ',I4.4,2X,'Ave= ',G11.5,1X,
+ 'Rms= ',G11.5,1X,'Ave2= ',G11.5,1X,'Rms2= ',G11.5)
ELSE
WRITE (*,40) IB,UBAR,URMS,VBAR,VRMS
40 FORMAT ('+','Rec ',I4.4,2X,'Ave= ',G11.5,1X,
+ 'Rms= ',G11.5,1X,'Ave2= ',G11.5,1X,'Rms2= ',G11.5)
ENDIF
* Save transformed values for this record.
IF (.NOT. NOFIL)
+ CALL FILTER (FTYPE,FC,F1,F2,NS,DELT,N,FR)
IF (RMEAN) THEN
DO I=1,N
FR(I)=FR(I)-UBAR
ENDDO
ENDIF
IF (DTREND) CALL DTRD (FR,N,DELT)
WRITE (NUMO,ERR=50) (FR(I), I=1,N)
IF (.NOT. NOFIL2)
+ CALL FILTER (FTYP2,FC2,F12,F22,NS2,DELT,N,FR2)
IF (RMEAN2) THEN
DO I=1,N
FR2(I)=FR2(I)-VBAR
ENDDO
ENDIF
IF (DTREND2) CALL DTRD (FR2,N,DELT)
WRITE (NUMO2,ERR=50) (FR2(I), I=1,N)
ENDDO
* All records of this data file have been processed
* and stored. Compute mean and rms values for this
* data file and continue.
TOTAL=FLOAT(NUMREC*N)
AVEU=SUM1UT/TOTAL
RMSU=SQRT(SUM2UT/TOTAL-AVEU**2)
AVEV=SUM1VT/TOTAL
RMSV=SQRT(SUM2VT/TOTAL-AVEV**2)
* Display the results for this file.
WRITE (*,'(/1X,A,2(A,G11.5,1X,A,G11.5,1X)/)')
+ 'File : ','Ave= ',AVEU,'Rms= ',RMSU,
+ 'Ave2= ',AVEV,'Rms2= ',RMSV
RETURN
50 WRITE (*,'(/13X,A,A/)') 'Run-time I/O error : ',
+ 'No such device or out of space !'
TRANS_TWO_CHAN=1
RETURN
END
FUNCTION WIN_SCALE (WIN)
IMPLICIT REAL*4 (A-H,O-Z), INTEGER*2 (I-N)
INCLUDE 'SPECTRUM.FN' ! Common statements
REAL*4 WIN[HUGE](*)
DATA A0,A1,A2,A3 /0.355768,0.487396,0.144232,0.012604/
DATA C0,C1,C2,C3 /0.3125,0.46875,0.1875,0.03125/
* Initialization.
WIN_SCALE=0
S=0.0
TNM1=FLOAT(N-1)
* Compute the scale factor.
DO I=1,N
TIM1=FLOAT(I-1)
IF (HANN)
+ WIN(I)=(1.0-COS(PI2*TIM1/TNM1))/2.0
IF (WELCH)
+ WIN(I)=1.0-((TIM1-TNM1/2.0)/(TNM1/2.0))**2
IF (PARZ)
+ WIN(I)=1.0-ABS((TIM1-TNM1/2.0)/(TNM1/2.0))
IF (TUKEY) THEN
IF (I .LE. K .OR. I .GE. N-K)
+ WIN(I)=(1.0-COS(10.0*PI*TIM1/TNM1))/2.0
IF (I .GT. K .AND. I .LT. N-K) WIN(I)=1.0
ENDIF
IF (MINSL)
+ WIN(I)=A0-A1*COS(PI2*TIM1/TNM1)
+ +A2*COS(2.0*PI2*TIM1/TNM1)-A3*COS(3.0*PI2*TIM1/TNM1)
IF (MXDEC)
+ WIN(I)=C0-C1*COS(PI2*TIM1/TNM1)
+ +C2*COS(2.0*PI2*TIM1/TNM1)-C3*COS(3.0*PI2*TIM1/TNM1)
S=S+WIN(I)*WIN(I)
ENDDO
S=S/FLOAT(N)
S=SQRT(1.0/S)
RETURN
40 WRITE (*,'(/13X,A,A/)') 'Run-time I/O error : ',
+ 'No such device or out of space !'
WIN_SCALE=1
RETURN
END
FUNCTION TAPER_DATA (FR,FR2,FI,FI2,WIN)
IMPLICIT REAL*4 (A-H,O-Z), INTEGER*2 (I-N)
INCLUDE 'SPECTRUM.FN' ! Common statements
REAL*4 FR[HUGE](*)
REAL*4 FR2[HUGE](*)
REAL*4 FI[HUGE](*)
REAL*4 FI2[HUGE](*)
REAL*4 WIN[HUGE](*)
* Initialization.
TAPER_DATA=0
IBO=0
* Do the tapering now.
DO IB=1,NUMREC+1
IF (IB .EQ. NUMREC+1 .AND. NOOVRLP) RETURN
* Reposition input file if last time thru loop.
IF (IB .EQ. NUMREC+1 .AND. OVRLP) THEN
REWIND (NUMO)
READ (NUMO) ICHANS,IRSIZE,NUMREC,IDELTMS
IF (TWOCHAN) THEN
REWIND (NUMO2)
READ (NUMO2) ICHANS,IRSIZE,NUMREC,IDELTMS
ENDIF
ENDIF
* Get the next N points.
READ (NUMO) (FR(I),I=1,N)
IF (TWOCHAN) READ (NUMO2) (FR2(I),I=1,N)
IF (IB .EQ. 1 .OR. NOOVRLP) GO TO 10
* Operate on the first half of FR and store in FI.
DO I=N2+1,N
FI(I)=FR(I-N2)*S*WIN(I)
IF (TWOCHAN) FI2(I)=FR2(I-N2)*S*WIN(I)
ENDDO
* Write contents of FI to output file.
IBO=IBO+1
WRITE (*,'(''+'',24X,''Tapering record '',I4.4)') IBO
WRITE (NUMT,ERR=20) (FI(I),I=1,N)
IF (TWOCHAN) WRITE (NUMT2,ERR=20) (FI2(I),I=1,N)
IF (IB .EQ. NUMREC+1) RETURN
* Operate on all of FR.
10 DO I=1,N
* Operate on the second half of FR and store in FI.
IF (I .LE. N2 .AND. OVRLP) THEN
FI(I)=FR(N2+I)*S*WIN(I)
IF (TWOCHAN) FI2(I)=FR2(N2+I)*S*WIN(I)
ENDIF
* Now operate on all of FR.
FR(I)=FR(I)*S*WIN(I)
IF (TWOCHAN) FR2(I)=FR2(I)*S*WIN(I)
ENDDO
* Write contents of FR to output file.
IBO=IBO+1
IF (IBO .EQ. 1) THEN
WRITE (*,'(25X,''Tapering record '',I4.4)') IBO
ELSE
WRITE (*,'(''+'',24X,''Tapering record '',I4.4)') IBO
ENDIF
WRITE (NUMT,ERR=20) (FR(I),I=1,N)
IF (TWOCHAN) WRITE (NUMT2,ERR=20) (FR2(I),I=1,N)
ENDDO
RETURN
20 WRITE (*,'(/13X,A,A/)') 'Run-time I/O error : ',
+ 'No such device or out of space !'
TAPER_DATA=1
RETURN
END
FUNCTION SPEC_ONE_CHAN (FR,FI,FR2,FI2)
IMPLICIT REAL*4 (A-H,O-Z), INTEGER*2 (I-N)
INCLUDE 'SPECTRUM.FN' ! Common statements
REAL*4 FR[HUGE](*)
REAL*4 FI[HUGE](*)
REAL*4 FR2[HUGE](*)
REAL*4 FI2[HUGE](*)
* Initialization.
SPEC_ONE_CHAN=0
* Spectrum computations for one channel.
DO IB=1,NRECD2
IF (IB .EQ. 1) THEN
WRITE (*,'(/25X,''Transforming records '',I4.4,
+ '' and '',I4.4)') 2*IB-1,2*IB
ELSE
WRITE (*,'(''+'',24X,''Transforming records '',I4.4,
+ '' and '',I4.4)') 2*IB-1,2*IB
ENDIF
IF ( (OVRLP .OR. ODDNUM) .AND. IB .EQ. NRECD2)
+ WRITE (*,'(/25X,A,I4.4,A)')
+ 'Record ',2*IB,' was discarded.'
* Process two records simultaneously. This increases speed but
* requires an even number of records in the one channel case.
READ (INUM) (FR(I),I=1,N)
IF (ODDNUM .AND. IB .EQ. NRECD2) THEN
DO I=1,N
FI(I)=FR(I)
ENDDO
ELSE
READ (INUM) (FI(I),I=1,N)
ENDIF
CALL FFT (FR,FI,IOD)
IF ( (OVRLP .OR. ODDNUM) .AND. IB .EQ. NRECD2) FI(1)=0.0
XMEAN=XMEAN+FR(1)
YMEAN=YMEAN+FI(1)
DO K=1,N2
AK=FR(K+1)
BK=FI(K+1)
ANK=FR(NP1-K)
BNK=FI(NP1-K)
* Unscramble X and Y.
XR=AK+ANK
XI=BK-BNK
YR=BNK+BK
YI=ANK-AK
IF (POWER) THEN
B1K=XR*XR+XI*XI
B2K=YR*YR+YI*YI
ELSE
B1K=SQRT(XR*XR+XI*XI)
B2K=SQRT(YR*YR+YI*YI)
ENDIF
IF ( (OVRLP .OR. ODDNUM) .AND. IB .EQ. NRECD2) B2K=0.0
FR2(K)=FR2(K)+B1K
FI2(K)=FI2(K)+B2K
ENDDO
ENDDO
CLOSE (INUM,STATUS='DELETE',ERR=10)
IB=IB-1 ! Decrement added count in loop variable
* Put results into one buffer.
XMEAN=XMEAN+YMEAN
DO I=1,N
FR2(I)=FR2(I)+FI2(I)
ENDDO
IB=IB*2
IF (OVRLP .OR. ODDNUM) IB=IB-1 ! Remove bad record.
* Finished.
FN=FLOAT(N)
TOTAL=FN*FLOAT(IB)
XMEAN=XMEAN/TOTAL
DT=DELT
DELF=1.0/(FN*DT)
IF (POWER) SCALE=DT/(2.0*TOTAL)
IF (AMP) THEN
SCALE=DT/(2.0*FLOAT(IB))
XMEAN=XMEAN*TOTAL*SCALE*2.0
ENDIF
DO I=1,N
FR2(I)=FR2(I)*SCALE
ENDDO
IF (POWER) XRMS=SQRT(DELF*ADDSUM(FR2,N2))
* Put results into output [.PRN] file.
IPD=INDEX (POUT,'.',.TRUE.)
WRITE (*,'(/25X,''Writing '',A,'' now.'')') POUT(:IPD+3)
OPEN (NUMP,FILE=POUT,STATUS='UNKNOWN',ERR=10)
IF (POWER) THEN
WRITE (NUMP,'(F11.5,2X,F11.5/)',ERR=10) XMEAN,XRMS
DO I=1,N2
IF (FR2(I) .LT. EPS) FR2(I)=EPS
WRITE (NUMP,'(F11.5,2X,G15.5,2X,G15.5,2X,G15.5)',ERR=10)
+ DELF*I,ALOG10(DELF*I),FR2(I),ALOG10(FR2(I))
ENDDO
ELSE
WRITE (NUMP,'(G15.5,2X,G15.5,2X,G15.5)',ERR=10)
+ EPS,ALOG10(EPS),XMEAN
WRITE (NUMP,'(F11.5,2X,G15.5,2X,G15.5)',ERR=10)
+ (DELF*I,ALOG10(DELF*I),FR2(I), I=1,N2)
ENDIF
CLOSE (NUMP,STATUS='KEEP',ERR=10)
* Display results on the screen.
IF (POWER) THEN
WRITE (*,'(/29X,A,10X,A/21X,1P,2G15.5)')
+ 'mean','rms',XMEAN,XRMS
PERR = 1.0/SQRT(FLOAT(NREC)) ! use numrec before overlap
IF (TAPER .AND. NOOVRLP) PERR=PERR*SQRT(2.0)
WRITE (*,'(/25X,A,F6.5)') 'norm. random error = ',PERR
ENDIF
RETURN
10 WRITE (*,'(/13X,A,A/)') 'Run-time I/O error : ',
+ 'No such device or out of space !'
SPEC_ONE_CHAN=1
RETURN
END
FUNCTION SPEC_TWO_CHAN (FR,FI,FR2,FI2,WIN)
IMPLICIT REAL*4 (A-H,O-Z), INTEGER*2 (I-N)
INCLUDE 'SPECTRUM.FN' ! Common statements
REAL*4 FR[HUGE](*)
REAL*4 FI[HUGE](*)
REAL*4 FR2[HUGE](*)
REAL*4 FI2[HUGE](*)
REAL*4 WIN[HUGE](*)
* Initialization.
SPEC_TWO_CHAN=0
* Spectrum computations for two channels.
DO IB=1,NUMREC
IF (IB .EQ. 1) THEN
WRITE (*,'(/25X,''Transforming record '',I4.4)') IB
ELSE
WRITE (*,'(''+'',24X,''Transforming record '',I4.4)') IB
ENDIF
IF (OVRLP .AND. IB .EQ. NUMREC)
+ WRITE (*,'(/25X,A,I4.4,A)')
+ 'Record ',IB+1, ' was discarded for each channel.'
READ (INUM) (FR(I),I=1,N)
READ (INUM2) (FI(I),I=1,N)
CALL FFT (FR,FI,IOD)
XMEAN=XMEAN+FR(1)
YMEAN=YMEAN+FI(1)
DO K=1,N2
AK=FR(K+1)
BK=FI(K+1)
ANK=FR(NP1-K)
BNK=FI(NP1-K)
* Unscramble X and Y.
XR=AK+ANK
XI=BK-BNK
YR=BNK+BK
YI=ANK-AK
IF (POWER) THEN
B1K=XR*XR+XI*XI
B2K=YR*YR+YI*YI
ELSE IF (CROSS) THEN
B1K=XR*YR+XI*YI
B2K=XI*YR-XR*YI
IF (COH) THEN
B3K=XR*XR+XI*XI
B4K=YR*YR+YI*YI
ENDIF
ELSE
B1K=SQRT(XR*XR+XI*XI)
B2K=SQRT(YR*YR+YI*YI)
ENDIF
FR2(K)=FR2(K)+B1K
FI2(K)=FI2(K)+B2K
IF (COH) THEN
WIN(K)=WIN(K)+B3K
WIN(N2+K)=WIN(N2+K)+B4K
ENDIF
ENDDO
ENDDO
CLOSE (INUM,STATUS='DELETE',ERR=20)
CLOSE (INUM2,STATUS='DELETE',ERR=20)
IB=IB-1 ! Decrement added count in loop variable
* Finished.
FN=FLOAT(N)
TOTAL=FN*FLOAT(IB)
XMEAN=XMEAN/TOTAL
YMEAN=YMEAN/TOTAL
DT=DELT
DELF=1.0/(FN*DT)
IF (POWER .OR. CROSS) SCALE=DT/(2.0*TOTAL)
IF (AMP) THEN
SCALE=DT/(2.0*FLOAT(IB))
XMEAN=XMEAN*TOTAL*SCALE*2.0
YMEAN=YMEAN*TOTAL*SCALE*2.0
ENDIF
DO I=1,N
FR2(I)=FR2(I)*SCALE
FI2(I)=FI2(I)*SCALE
ENDDO
IF (COH) THEN
DO I=1,N
WIN(I)=WIN(I)*SCALE
ENDDO
ENDIF
IF (POWER) THEN
XRMS=SQRT(DELF*ADDSUM(FR2,N2))
YRMS=SQRT(DELF*ADDSUM(FI2,N2))
ENDIF
IF (CROSS) THEN
DEGRAD=360.0/PI2
S2=SQRT(2.0)
NR=NREC ! use numrec before overlap
SNR=SQRT(FLOAT(NR))
S2NR=SQRT(FLOAT(NR)*2.0)
RXY0=DELF*ADDSUM(FR2,N2)
DO I=1,N2
FII=FLOAT(I)
GXYMAG=SQRT(FR2(I)**2+FI2(I)**2)
PHASE=ATAN2(FI2(I),FR2(I))
FR2(I)=GXYMAG
FI2(I)=PHASE
FI2(N2+I)=PHASE
* If coherence is calculated then compute errors associated
* with cross spectrum estimates.
IF (COH) THEN
GAMMA=FR2(I)**2/(WIN(I)*WIN(N2+I))
WIN(I)=GAMMA
FR(I)=1.0/(SQRT(WIN(I))*SNR) ! e [ |Gxy(f)| ]
FR(N2+I)=SQRT(1.0-WIN(I))/(SQRT(WIN(I))*S2NR) ! s.d. [ Oxy(f) ]
FR(N2+I)=FR(N2+I)*DEGRAD ! s.d. [O] in deg
FI(I)=S2*(1.0-WIN(I))/(SQRT(WIN(I))*SNR) ! e [ gxy2(f) ]
IF (TAPER .AND. NOOVRLP) THEN
FR(I)=FR(I)*S2
FR(N2+I)=FR(N2+I)*S2
FI(I)=FI(I)*S2
ENDIF
ENDIF
ENDDO
* Apply a phase shift to straighten out Theta only if
* coherence is calculated -- need s.d. estimate.
IF (COH) THEN
NSD2=(NSD-1)/2
* Check on proper relationship of NPTS and NSD.
IF (NSD2+1 .GT. NPTS) THEN
WRITE (*,'(/21X,A,A/21X,A)') 'NSD too big or NPTS',
+ ' too small.','No adjustment of phase angle performed.'
GO TO 10
ENDIF
* Get first NPTS points. Pass this to fit subroutine.
DO I=1,NPTS
X(I)=DELF*I
Y(I)=FI2(N2+I)
ENDDO
* Sum first NSD points to check on average s.d. of theta.
SDSUM=0.0
DO I=1,NSD
SDSUM=SDSUM+FR(N2+NPTS-NSD2-1+I)
ENDDO
DO I=NPTS+1,N2
* Check on moving average of s.d. of theta.
* If too high -- quit.
SDSUM=SDSUM+FR(N2+I+NSD2)
SDSUM=SDSUM-FR(N2+I-NSD2-1)
SDAVG=SDSUM/NSD
IF (SDAVG .GT. SDLIM) GO TO 10
* Perform linear least squares fit on NPTS points before
* the current point. Guess value of current point from fit.
CALL LSF (X,Y,NPTS,XA1,XB2)
TEST=XA1+XB2*DELF*I
* XLO & XHI are current point guess +/- pi.
XLO=TEST-PI
XHI=TEST+PI
* If actual value is .GT. pi away from guessed value then
* add 2*pi to or subtract 2*pi from actual value as required.
DO WHILE (FI2(N2+I) .LT. XLO)
FI2(N2+I)=FI2(N2+I)+PI2
ENDDO
DO WHILE (FI2(N2+I) .GT. XHI)
FI2(N2+I)=FI2(N2+I)-PI2
ENDDO
* Prepare set of NPTS points before next current point.
DO J=1,NPTS
X(J)=DELF*(I-NPTS+J)
Y(J)=FI2(N2+I-NPTS+J)
ENDDO
ENDDO
ENDIF
ENDIF
* Put results into output [.PRN] file.
10 IPD=INDEX (POUT,'.',.TRUE.)
WRITE (*,'(/25X,''Writing '',A,'' now.''/)') POUT(:IPD+3)
OPEN (NUMP,FILE=POUT,STATUS='UNKNOWN',ERR=20)
IF (POWER) THEN
WRITE (NUMP,'(4(F11.5,2X))',ERR=20) XMEAN,XRMS,YMEAN,YRMS
WRITE (NUMP,'( )',ERR=20)
DO I=1,N2
IF (FR2(I) .LT. EPS) FR2(I)=EPS
IF (FI2(I) .LT. EPS) FI2(I)=EPS
WRITE (NUMP,'(F11.5,1X,F11.5,1X,G15.5,1X,
+ F11.5,1X,G15.5,1X,F11.5)',ERR=20)
+ DELF*I,ALOG10(DELF*I),FR2(I),ALOG10(FR2(I)),
+ FI2(I),ALOG10(FI2(I))
ENDDO
ELSE IF (CROSS) THEN
IF (NOCOH) THEN
WRITE (NUMP,'(G15.5)',ERR=20) RXY0
WRITE (NUMP,'( )',ERR=20)
WRITE (NUMP,'(F11.5,1X,F11.5,1X,G15.5,1X,F11.5,
+ 1X,F11.5)',ERR=20) (DELF*I,ALOG10(DELF*I),
+ FR2(I),ALOG10(FR2(I)),FI2(I),I=1,N2)
ELSE IF (COH) THEN
WRITE (NUMP,'(G15.5)',ERR=20) RXY0
WRITE (NUMP,'( )',ERR=20)
WRITE (NUMP,'(F11.5,1X,F11.5,1X,G15.5,1X,F11.5,1X,
+ F11.5,1X,F11.5,1X,G15.5)',ERR=20) (DELF*I,
+ ALOG10(DELF*I),FR2(I),ALOG10(FR2(I)),FI2(I),FI2(N2+I),
+ WIN(I),I=1,N2)
ENDIF
ELSE
WRITE (NUMP,'(4(G15.5,2X))',ERR=20)
+ EPS,ALOG10(EPS),XMEAN,YMEAN
WRITE (NUMP,'(F11.5,2X,G15.5,2X,G15.5,2X,G15.5)',ERR=20)
+ (DELF*I,ALOG10(DELF*I),FR2(I),FI2(I),I=1,N2)
ENDIF
CLOSE (NUMP,STATUS='KEEP',ERR=20)
* Display results on the screen.
IF (POWER) THEN
WRITE (*,'(29X,A,10X,A,2(/21X,1P,2G15.5))')
+ 'mean','rms',XMEAN,XRMS,YMEAN,YRMS
PERR = 1.0/SQRT(FLOAT(NREC)) ! use numrec before overlap
IF (TAPER .AND. NOOVRLP) PERR=PERR*SQRT(2.0)
WRITE (*,'(/25X,A,F6.5)') 'norm. random error = ',PERR
ENDIF
IF (CROSS) WRITE (*,'(21X,A,1P,G15.5)')
+ 'Rxy(0) = E [x(t) y(t)] = ',RXY0
* If coherence was calculated then store
* error estimates into separate error file.
IF (COH) THEN
IPD=INDEX (EOUT,'.',.TRUE.)
WRITE (*,'(/25X,''Writing '',A,'' now.'')') EOUT(:IPD+3)
OPEN (NUME,FILE=EOUT,STATUS='UNKNOWN',ERR=20)
WRITE (NUME,'(5(F11.5,2X))',ERR=20) (DELF*I,
+ ALOG10(DELF*I),FR(I),FR(N2+I),FI(I),I=1,N2)
CLOSE (NUME,STATUS='KEEP',ERR=20)
ENDIF
RETURN
20 WRITE (*,'(/13X,A,A/)') 'Run-time I/O error : ',
+ 'No such device or out of space !'
SPEC_TWO_CHAN=1
RETURN
END
********************************************************
SUBROUTINE FFT (FR,FI,K)
* *
* Fast Fourier Transform of a *
* Complex Time Series *
* Features : time decomposition (odd-even) *
* : input bit reversal *
* : computation 'in place' *
* : N= number of points *
* : FR+j*FI - complex series *
* Source : Stearns 1975 *
********************************************************
IMPLICIT INTEGER*4 (I-N)
REAL*4 FR[HUGE](*)
REAL*4 FI[HUGE](*)
N=2**K
MR=0
NN=N-1
* Input bit reversal.
DO M=1,NN
L=N
10 L=L/2
IF(MR+L .GT. NN) GO TO 10
MR=MOD(MR,L)+L
IF (MR.LE. M) CYCLE
TR=FR(M+1)
FR(M+1)=FR(MR+1)
FR(MR+1)=TR
TI=FI(M+1)
FI(M+1)=FI(MR+1)
FI(MR+1)=TI
ENDDO
* Compute the coefficients.
L=1
PI=2.0*ASIN(1.0)
20 IF (L.GE.N) RETURN
ISTEP=2*L
EL=L
DO M=1,L
A=PI*FLOAT(1-M)/EL
WR=COS(A)
WI=SIN(A)
DO I=M,N,ISTEP
J=I+L
TR=WR*FR(J)-WI*FI(J)
TI=WR*FI(J)+WI*FR(J)
FR(J)=FR(I)-TR
FI(J)=FI(I)-TI
FR(I)=FR(I)+TR
FI(I)=FI(I)+TI
ENDDO
ENDDO
L=ISTEP
GO TO 20
END
SUBROUTINE RDCFG (ROOT1,ROOT2,SUFX,PTH,TPTH,SUFFIX)
* Routine to read config file - spectrum.cfg
PARAMETER (NUM=2)
LOGICAL*1 ROOT1,ROOT2,SUFX
CHARACTER*1 FIRST
CHARACTER*4 SUFFIX
CHARACTER*72 FLINE,PTH,TPTH
* File exists. Get various paths.
KOUNT=0
OPEN (NUM,FILE='SPECTRUM.CFG',STATUS='OLD')
IF (EOF(NUM)) GO TO 30
10 READ (NUM,'(A)',ERR=30) FLINE
* Check for a comment line.
IF (FLINE (1:1) .EQ. '*') GO TO 10
KOUNT=KOUNT+1
* Convert to upper case characters if necessary.
DO J=72,1,-1
FIRST=FLINE(J:J)
IF (ICHAR(FIRST) .GE. 97 .AND. ICHAR(FIRST) .LE. 122) THEN
IHOLD=ICHAR(FIRST)-32
FIRST=CHAR(IHOLD)
FLINE(J:J)=FIRST
ENDIF
ENDDO
* Assign appropriate paths.
IF (KOUNT .EQ. 1) THEN
IF (FLINE(:4) .EQ. 'ROOT') THEN
ROOT1=.TRUE.
ELSE
ROOT1=.FALSE.
ISL=INDEX (FLINE,'\',.TRUE.)
PTH=FLINE(:ISL)
ENDIF
ELSE IF (KOUNT .EQ. 2) THEN
IF (FLINE(:4) .EQ. 'ROOT') THEN
ROOT2=.TRUE.
ELSE
ROOT2=.FALSE.
ISL=INDEX (FLINE,'\',.TRUE.)
TPTH=FLINE(:ISL)
ENDIF
ELSE IF (KOUNT .EQ. 3) THEN
SUFX=.TRUE.
SUFFIX=FLINE(:4)
ENDIF
GO TO 10
30 CLOSE (NUM)
RETURN
END
SUBROUTINE LSF(X,Y,N,A,B)
IMPLICIT REAL*4 (A-H,O-Z), INTEGER*2 (I-N)
REAL*4 X(N),Y(N)
* Initialization.
XS=0.0
X2S=0.0
YS=0.0
XYS=0.0
* Compute sums.
DO I=1,N
XS=XS+X(I)
X2S=X2S+X(I)*X(I)
YS=YS+Y(I)
XYS=XYS+X(I)*Y(I)
ENDDO
* Compute coefficients and guess next value.
DEN=N*X2S-XS*XS
A=(X2S*YS-XS*XYS)/DEN
B=(N*XYS-XS*YS)/DEN
RETURN
END
FUNCTION ADDSUM (Y,N)
* SUMMING of elements of real*4 array Y. This function uses
* a tree summming process to minimize truncation error for
* large N(>500). A Real*4 and Integer*2 buffer of size log2(n)
* is required. This version received from Dr. Anuvat Sirivat.
REAL*4 Y[HUGE](*)
REAL*4 B(16),ADDSUM
INTEGER*2 W(16)
INTEGER*4 N
K=1
W(1)=0
B(1)=0.0
DO I=1,N
B(K)=B(K)+Y(I)
W(K)=W(K)+1
10 K1=K-1
IF(K1.LE.0) GO TO 20
IF(W(K1).GT.W(K)) GO TO 20
B(K1)=B(K1)+B(K)
W(K1)=W(K1)+W(K)
K=K-1
GO TO 10
20 IF(W(K).LT.2) CYCLE
K=K+1
W(K)=0
B(K)=0.0
ENDDO
30 K2=K/2
IF(K2.EQ.0) GO TO 40
J=1
DO I=1,K2
B(I)=B(J)+B(J+1)
J=J+2
ENDDO
IF(J.EQ.K) B(1)=B(1)+B(K)
K=K2
GO TO 30
40 ADDSUM=B(1)
RETURN
END
SUBROUTINE DTRD (FR,N,DELT)
* This routine will remove a linear mean trend from the
* original data by means of least squares fit. This is a
* modified version of a routine received from Dr. Sirivat.
IMPLICIT REAL*4 (A-H,O-Z), INTEGER*4 (I-N)
INTEGER*4 N
REAL*4 DELT
REAL*4 FR[HUGE](*)
* Initialization.
FN=FLOAT(N)
YSR=0.0
XYSR=0.0
* Compute required sums for least squares fit.
DO I=1,N
FII=FLOAT(I)
YSR=YSR+FR(I)
XYSR=XYSR+FII*FR(I)
ENDDO
XS=FN*(FN+1.0)/2.0*DELT
X2S=FN*(FN+1.0)*(2.0*FN+1.0)/6.0*DELT*DELT
XYSR=XYSR*DELT
DER=FN*X2S-XS*XS
* Compute coefficients and remove trend.
B0R=(YSR*X2S-XYSR*XS)/DER
B1R=(FN*XYSR-XS*YSR)/DER
DO I=1,N
FII=FLOAT(I)
FR(I)=FR(I)-B0R-B1R*FII*DELT
ENDDO
RETURN
END
SUBROUTINE FILTER (FTYPE,FC,F1,F2,NS,DELT,N,FR)
* Routine to filter a data file. (FILTER.FOR)
* This routine will filter a data file using any one of
* four different types of IIR recursive digital filters:
* Lowpass, Highpass, Bandpass, and Bandstop. The routines
* are modified versions of those given in Stearns, 1975.
PARAMETER (NSMAX=20)
LOGICAL*1 LPASS,HPASS,BPASS,BSTOP
INTEGER*2 NS
INTEGER*4 N,J
REAL*4 FR[HUGE](*)
REAL*4 X(NSMAX+1,5)
REAL*4 A(NSMAX),B(NSMAX),C(NSMAX),D(NSMAX),E(NSMAX),G
CHARACTER*1 FTYPE
* Initialization.
X=0.0
G=0.0
LPASS=(FTYPE .EQ. 'L')
HPASS=(FTYPE .EQ. 'H')
BPASS=(FTYPE .EQ. 'B')
BSTOP=(FTYPE .EQ. 'S')
* Design the desired filter.
IF (LPASS .OR. HPASS) KHI=3
IF (BPASS .OR. BSTOP) KHI=5
IF (LPASS) CALL LP (FC,DELT,NS,A,B,C)
IF (HPASS) CALL HP (FC,DELT,NS,A,B,C)
IF (BPASS) CALL BP (F1,F2,DELT,NS,A,B,C,D,E)
IF (BSTOP) CALL BS (F1,F2,DELT,NS,A,B,C,D,E,G)
G2=2.0*G
GG=2.0+G*G
* Filter the array FR(J).
DO J=1,N
* Go thru NS sections of filter.
IF (LPASS) THEN
X(1,3)=FR(J)
DO I=1,NS
TEMP=A(I)*(X(I,3)+2.0*X(I,2)+X(I,1))-B(I)*X(I+1,2)
X(I+1,3)=TEMP-C(I)*X(I+1,1)
ENDDO
ELSE IF (HPASS) THEN
X(1,3)=FR(J)
DO I=1,NS
TEMP=A(I)*(X(I,3)-2.0*X(I,2)+X(I,1))-B(I)*X(I+1,2)
X(I+1,3)=TEMP-C(I)*X(I+1,1)
ENDDO
ELSE IF (BPASS) THEN
X(1,5)=FR(J)
DO I=1,NS
TEMP=A(I)*(X(I,5)-2.0*X(I,3)+X(I,1))-B(I)*X(I+1,4)
X(I+1,5)=TEMP-C(I)*X(I+1,3)-D(I)*X(I+1,2)-E(I)*X(I+1,1)
ENDDO
ELSE IF (BSTOP) THEN
X(1,5)=FR(J)
DO I=1,NS
TEMP=A(I)*(X(I,5)+G2*X(I,4)+GG*X(I,3)+G2*X(I,2)+X(I,1))
X(I+1,5)=TEMP-B(I)*X(I+1,4)-C(I)*X(I+1,3)-D(I)*X(I+1,2)
+ -E(I)*X(I+1,1)
ENDDO
ENDIF
* Push down past values in filter.
DO I=1,NS+1
DO K=1,KHI-1
X(I,K)=X(I,K+1)
ENDDO
ENDDO
* Set FR(J) equal to the output of the filter and continue.
FR(J)=X(NS+1,KHI)
ENDDO
RETURN
END
SUBROUTINE LP (FC,T,NS,A,B,C)
* Lowpass Butterworth Digital Filter
* Design Subroutine (LP.FOR)
INTEGER*2 NS
REAL*4 A(*),B(*),C(*)
PI=2.0*ASIN(1.0)
WCP=SIN(FC*PI*T)/COS(FC*PI*T)
DO K=1,NS
CS=COS(FLOAT(2*(K+NS)-1)*PI/FLOAT(4*NS))
X=1.0/(1.0+WCP*WCP-2.0*WCP*CS)
A(K)=WCP*WCP*X
B(K)=2.0*(WCP*WCP-1.0)*X
C(K)=(1.0+WCP*WCP+2.0*WCP*CS)*X
ENDDO
RETURN
END
SUBROUTINE HP (FC,T,NS,A,B,C)
* Highpass Butterworth Digital Filter
* Design Subroutine (HP.FOR)
INTEGER*2 NS
REAL*4 A(*),B(*),C(*)
PI=2.0*ASIN(1.0)
WCP=SIN(FC*PI*T)/COS(FC*PI*T)
DO K=1,NS
CS=COS(FLOAT(2*(K+NS)-1)*PI/FLOAT(4*NS))
A(K)=1.0/(1.0+WCP*WCP-2.0*WCP*CS)
B(K)=2.0*(WCP*WCP-1.0)*A(K)
C(K)=(1.0+WCP*WCP+2.0*WCP*CS)*A(K)
ENDDO
RETURN
END
SUBROUTINE BP (F1,F2,T,NS,A,B,C,D,E)
* Bandpass Butterworth Digital Filter
* Design Subroutine (BP.FOR)
INTEGER*2 NS
REAL*4 A(*),B(*),C(*),D(*),E(*)
PI=2.0*ASIN(1.0)
W1=SIN(F1*PI*T)/COS(F1*PI*T)
W2=SIN(F2*PI*T)/COS(F2*PI*T)
WC=W2-W1
Q=WC*WC+2.0*W1*W2
S=W1*W1*W2*W2
DO K=1,NS
CS=COS(FLOAT(2*(K+NS)-1)*PI/FLOAT(4*NS))
P=-2.0*WC*CS
R=P*W1*W2
X=1.0+P+Q+R+S
A(K)=WC*WC/X
B(K)=(-4.0-2.0*P+2.0*R+4.0*S)/X
C(K)=(6.0-2.0*Q+6.0*S)/X
D(K)=(-4.0+2.0*P-2.0*R+4.0*S)/X
E(K)=(1.0-P+Q-R+S)/X
ENDDO
RETURN
END
SUBROUTINE BS (F1,F2,T,NS,A,B,C,D,E,G)
* Bandstop Butterworth Digital Filter
* Design Subroutine (BS.FOR)
INTEGER*2 NS
REAL*4 A(*),B(*),C(*),D(*),E(*),G
PI=2.0*ASIN(1.0)
W1=SIN(F1*PI*T)/COS(F1*PI*T)
W2=SIN(F2*PI*T)/COS(F2*PI*T)
WC=W2-W1
XK=1.0+W1*W2
G=(2.0*XK-4.0)/XK
Q=WC*WC+2.0*W1*W2
S=W1*W1*W2*W2
DO K=1,NS
CS=COS(FLOAT(2*(K+NS)-1)*PI/FLOAT(4*NS))
P=-2.0*WC*CS
R=P*W1*W2
X=1.0+P+Q+R+S
A(K)=XK*XK/X
B(K)=(-4.0-2.0*P+2.0*R+4.0*S)/X
C(K)=(6.0-2.0*Q+6.0*S)/X
D(K)=(-4.0+2.0*P-2.0*R+4.0*S)/X
E(K)=(1.0-P+Q-R+S)/X
ENDDO
RETURN
END